function new_img=proj_distortion_removal(filename)
img = my_image_read(filename,2);
width=size(img,1);
height=size(img,2);

dim = [width height 1]';
%get input
figure(1);imshow(img)
for i=1:4
    [p(1,i) p(2,i)]=ginput(1);
    if (i>1)
            line = ones(size(img(:,:,1)));
            line = bresenham(line,[p(1,i-1) p(2,i-1);p(1,i) p(2,i)]);
            line = logical(1-line);
            for j=1:3
                layer = img(:,:,j);
                layer(line)=0;
                img(:,:,j)=layer;
            end
            figure(1);
            imshow(img);
    end
    
end
p(3,:)=1;

temp = p(:,3);
p(:,3) = p(:,4);
p(:,4) = temp;
    

l1=cross(p(:,1),p(:,2));
l2=cross(p(:,3),p(:,4));
m1=cross(p(:,1),p(:,3));
m2=cross(p(:,2),p(:,4));

v1 = cross(l1,l2);
v2 = cross(m1,m2);

v_line = cross(v1,v2);

v_line = v_line./v_line(end);
v_line(1:2)=v_line(1:2);
v_line(end)=v_line(end);

 H1 = [ 1 0 0;
        0 1 0;
        v_line'];   

img2 = myWarp(img,H1);
figure(2);
 imshow(img2);
 
 %%%%%%%%%%%%%%
 % Remove affine distortion
 %%%%%%%%%%%%%%
 
 %apply H1 on p
 
 p_n=H1*p;
 for i=1:4
     p_n(:,i)=p_n(:,i)/p_n(end,i);
 end
 l1=cross(p_n(:,1),p_n(:,2));
l2=cross(p_n(:,1),p_n(:,4));
m1=cross(p_n(:,1),p_n(:,3));
m2=cross(p_n(:,3),p_n(:,2));

M(1,:)=[l1(1)*m1(1) l1(1)*m1(2)+l1(2)*m1(1)];
M(2,:)=[l2(1)*m2(1) l2(1)*m2(2)+l2(2)*m2(1)];

b=-[l1(2)*m1(2) l2(2)*m2(2)]';
%solve M*x=b:
x=inv(M)*b;
%S=A*A'
S=[x(1) x(2);x(2) 1];
[u s v]=svd(S);
A=u*sqrt(s)*u';

H2=[A(1,1) A(1,2) 0;
        A(2,1) A(2,2) 0;
        0 0 1];
   
img3=myWarp(img,inv(H2)*H1);
figure (3);
imshow(img3);
new_img=img3;